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Precise understanding of strongly interacting fermions, from electrons in modern materials to 
nuclear matter, presents a major goal in modern physics. However, the theoretical description of 
interacting Fermi systems is usually plagued by the intricate quantum statistics at play. Here we 
present a cross-validation between a new theoretical approach. Bold Diagrammatic Monte Carlo 
(BDMC), and precision experiments on ultra-cold atoms. Specifically, we compute and measure 
with unprecedented accuracy the normal-state equation of state of the unitary gas, a prototypical 
example of a strongly correlated fermionic system. Excellent agreement demonstrates that a series 
of Feynman diagrams can be controUably resummed in a non-perturbative regime using BDMC. 
This opens the door to the solution of some of the most challenging problems across many areas of 
physics. 



In his seminal 1981 lecture [T], Feynman argued that 
an arbitrary quantum system cannot be efficiently simu- 
lated with a classical universal computer, because gener- 
ally, quantum statistics can only be imitated with a clas- 
sical theory if probabilities are replaced with negative (or 
complex) weighting factors. For the majority of many- 
particle models this indeed leads to the so-called sign 
problem which has remained an insurmountable obsta- 
cle. According to Feynman, the only way out is to employ 
computers made out of quantum-mechanical elements |T] . 
The recent experimental breakthroughs in cooling, prob- 
ing and controlling strongly interacting quantum gases 
prompted a challenging effort to use this new form of 
quantum matter to realize Feynman's emulators of fun- 
damental microscopic models [HH]. Somewhat ironically, 
Feynman's arguments which led him to the idea of em- 
ulators may be defied by a theoretical method that he 
himself devised, namely Feynman diagrams. This tech- 
nique organizes the calculation of a given physical quan- 
tity as a series of diagrams representing all possible ways 
particles can propagate and interact (see, e.g., Ref. [3]). 
For the many-body problem, this diagrammatic expan- 
sion is commonly used either in perturbative regimes or 
within uncontrolled approximations. However, the in- 
troduction of Diagrammatic Monte Carlo recently al- 
lowed to go well beyond the first few diagrams, and even 
reach convergence of the series in a moderately correlated 
regime [11 [5]. 

In this Letter we show that for a strongly correlated 
system and down to a phase transition, the diagrammatic 
series can still be given a mathematical meaning and 
leads to controllable results within Bold Diagrammatic 



Monte Carlo. This approach, proposed in Refs. [4l|6l[2], 
is first implemented here for the many-body problem. 
We focus on the unitary gas, i.e. spin-1/2 fermions with 
zero-range interactions at infinite scattering length [S]- 
flUj . This system offers the unique possibility to strin- 
gently test our theory against a quantum emulator real- 
ized here with trapped ultracold ^Li atoms at a broad 
Feshbach resonance [51-fTU] . This experimental validation 
is indispensable for our theory based on resummation of a 
possibly divergent series: although the physical answer is 
shown to be independent from the applied resummation 
technique ~ suggesting that the procedure is adequate 
- its mathematical validity remains to be proven. In 
essence. Nature provides the 'proof. This presents the 
first - though long-anticipated - compelling example of 
how ultra-cold atoms can guide new microscopic theories 
for strongly interacting quantum matter. 

At unitarity, the disappearance of an interaction- 
imposed length scale leads to scale invariance. This prop- 
erty renders the model relevant for other physical systems 
such as neutron matter. It also makes the balanced (i.e., 
spin-unpolarised) unitary gas ideally suited for the ex- 
perimental high-precision determination of the equation 
of state (EOS) described below. Finally, it implies the 
absence of a small parameter, making the problem noto- 
riously difficult to solve. 

In traditional Monte Carlo approaches, which simulate 
a finite piece of matter, the sign problem causes an expo- 
nential increase of the computing time with system size 
and inverse temperature. In contrast, BDMC simulates a 
mathematical answer in the thermodynamic limit. This 
radically changes the role of the fermionic sign. Dia- 
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FIG. 1. Bold Diagrammatic Monte Carlo evaluates skele- 
ton diagrammatic series for the self-energy E and the pair 
self-energy 11 (lower box). The diagrams are built on dressed 
one-body propagators G and pair propagators F, which them- 
selves are the solution of the Dyson and Bethe-Salpeter equa- 
tions (upper box). This cycle is repeated until convergence is 
reached. is the non-interacting propagator, and F" is the 
partially dressed pair propagator obtained by summing the 
bare ladder diagrams. 



grammatic contributions are sign-alternating with order, 
topology and values of internal variables. Since the num- 
ber of graphs grows factorially with diagram order, a 
near cancellation between these contributions is actually 
necessary for the series to be resummable by techniques 
requiring a finite radius of convergence. We find that this 
'sign blessing' indeed takes place. 

In essence, BDMC solves the full quantum many-body 
problem by stochastically summing all the skeleton di- 
agrams for irreducible single-particle self-energy S and 
pair self-energy 11, expressed in terms of bold (i.e., fully 
dressed) single-particle and pair propagators G and F 
which are determined self-consistently (see Fig. 1). The 
density EOS (i.e., the relation between total density n, 
chemical potential and temperature T) is given by G 
at zero distance and imaginary time, n(fi,T) — 2G{r = 
0, r = 0~). The thermodynamic limit can be taken ana- 
lytically. The sum of ladder diagrams built on the bare 
single-particle propagator defines a partially dressed pair 
propagator F". Since F" is well defined for the zero- 
range continuous-space interaction, the zero-range limit 
can also be taken analytically. This is in sharp contrast 
with other numerical methods pTHT5] where taking the 
thermodynamic and zero-range limits is computationally 
very expensive. BDMC performs a random walk in the 
space of irreducible diagrams using local updates. The 
simulation is run in a self-consistent cycle (along the lines 
of Ref. l&i) until convergence is reached. Full details will 
be presented elsewhere. In essence, our approach up- 
grades the standard many-body theories based on one 



lowest-order diagram (see, e.g., Refs. [HIITS] ) to millions 
of graphs. 

In the quantum degenerate regime, we do not observe 
convergence of the diagrammatic series for T, and 11 eval- 
uated up to order 9. Here, order N means E-diagrams 
with N vertices (i.e., N F-lines) and H-diagrams with 
iV — 1 vertices. To extract the infinite-order result, we 
apply the following Abelian resummation methods [16j . 
The contribution of all diagrams of order N is multi- 
plied by e"*^^""^ where A„ depends on the resummation 
method: (i) A„ — n logn (with Ao = 0) for Lindelof TB], 
(ii) A„ — {n—1) log(n— 1) (with Aq = Ai = 0) for "shifted 
Lindelof" , or (iii) A„ = n? for Gaussian [17] . A full sim- 
ulation is performed for each e, and the final result is 
obtained by extrapolating to e = 0, see Fig. 2. 

This protocol relies on the following crucial mathemat- 
ical assumptions; (i) the iV-th order contribution of the 
diagrammatic expansion for S (for fixed external vari- 
ables) is the iV-th coefficient of the Taylor series at z = 
of a function g{z) which has a non-zero convergence ra- 
dius, (ii) the analytic continuation g{\), performed by 
the above resummation methods |161 117) . is the physi- 
cally correct value of S. The same assumptions should 
hold for n. 

Proving these assumptions is an open mathematical 
challenge. Note that Dyson's collapse argument [TH] 
is not applicable to immediately disprove the assump- 
tion (i) oi a, non-zero convergence radius: Indeed, unlike 
QED, our skeleton series is not an expansion in powers 
of a coupling constant whose sign change would lead to 
an instability. A first important evidence for the validity 
of our mathematical assumptions is that the three differ- 
ent resummation methods yield consistent results. For 
an independent test, we turn to experiments. 

The present experiment furnishes high-precision data 
for the density n as a function of the local value V of the 
trapping potential (see Fig. 3 and Methods). We start 
the process by obtaining the EOS at high temperatures 
in the non-degenerate wings of the atom cloud where the 
virial expansion is applicable. Once the temperature and 
the chemical potential have been determined from fits 
to the wings of the cloud, the data closer to the cloud 
center provides a new prediction of the EOS. The process 
is iterated to access lower temperatures. 

Scale invariance allows to write the density EOS as 
n{fi,T)X^ = fiPn), with A = / [mkBT) the ther- 

mal de Broglie wavelength, /? = l/{kBT) the inverse tem- 
perature and / a universal function. A convenient nor- 
malization of the data is provided by the EOS of a non- 
interacting Fermi gas, uqX^ — /o(/3/^). In Fig. 4a we thus 
report the ratio n(/i,r)/no(M,r) = /(/3^)//o(/3^), bring- 
ing out the difference between the ideal and the strongly 
interacting Fermi gas. The Gibbs-Duhem relation allows 
us to also calculate the pressure at a given chemical po- 
tential, P{^^o, T) = J!^l dfx n{fx, T) = ^F{Pno), where 
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FIG. 2. Cross-validation between resummation proce- 
dure and experiment at — +1. Bold Diagrammatic 
Monte Carlo data for the dimensionless density nX^, as a 
function of the parameter e controlling the resummation pro- 
cedure, for three different resummation methods: Lindelof 
(blue circles), shifted Lindelof (black diamonds), and Gauss 
(open green squares). The solid lines are linear fits to the 
Monte Carlo data, their e — >■ extrapolation agrees within er- 
ror bars with the experimental data point (filled red square). 
[In the opposite limit e — >■ oo, the Lindelof (resp. shifted Lin- 
delof) curves will asymptote to the first [151 119] (resp. third) 
order results, shown by the dashed (resp. dash-dotted) line.] 
Error bars for each e represent the statistical error, together 
with the estimated systematic error coming from not sampling 
diagrams of order > 9. 
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FIG. 3. Constructing the EOS from in situ imaging. 
The atom cloud shown contains N = 8 x 10"* atoms for each 
spin-state, with a local Fermi energy of Ef ~ 370 nK at the 
center, a) Absorption image of the atomic cloud after quad- 
rant averaging, b) Reconstructed local density n{p,z). c) 
Equipotential averaging produces a low-noise density profile 
n vs V. Thermometry is performed by fitting the experi- 
mental data (red) to the known portion of the EOS (solid 
blue line), starting with the virial expansion for j3fj. < —1.25 
(green dashed line). In this example, the EOS is known 
to Pfi < —0.25, and the fit to the density profile yields 
T = 113 nK, and — 1.63. d) Given ^ and T, the den- 
sity profile can be rescaled to produce the EOS nX^ vs P/i. 



F{x) — dx'f{x'). We normalize it by the pressure of 
the ideal Fermi gas and show F(/3^)/Fo(/3/x), see Fig. 4b. 
The agreement between BDMC and experiment is ex- 
cellent. The comparison is sufhciently sensitive to vali- 
date the procedure of resumming and extrapolating (see 
Fig. 2). The result was checked to be independent of the 
maximal sampled diagram order iVmax G {7; 8; 9} within 
the error bars displayed in Fig. 2 for each e. The BDMC 
final error bar in Fig. 4 is the sum of the conservatively 
estimated systematic errors from the uncertainty of the 
e — ?> extrapolation and from the dependence on numeri- 
cal grids and cutoffs, the latter being reduced by analyti- 
cally treating high-momentum short-time singular parts. 
The systematic error in the experiment is determined to 
be about 1% by the independent determination of the 
EOS of the non-interacting Fermi gas. The experimen- 
tal error bars of Fig. 4 also include the statistical error, 
which is < 0.5% thanks to the scale invariance of the 
balanced unitary gas: irrespective of shot-to-shot fluctu- 
ations of atom number and temperature, all experimen- 
tal profiles contribute to the same scaled EOS-function 
/. The dominant uncertainty on the experimental EOS 
stems from the uncertainty in the position of the ^Li Fes- 
hbach resonance known to be at 834.15±1.5 G from spec- 
troscopic measurements [20j . The change in energy, pres- 
sure and density with respect to the interaction strength 
is controlled by the so-called contact [21] that is obtained 



from r in the BDMC calculation. This allows us to de- 
fine the uncertainty margins above and below the exper- 
imental data (see Fig. 4) that give the prediction for the 
unitary EOS if the true Feshbach resonance lied 1.5 G 
below or above 834.15 G, respectively. 

We clearly discriminate against previous theoretical 
and experimental results. Deviations from the theory 
based on the first-order Feynman diagrams [TSJ [TH] are 
expected, and rather remarkably moderate. Differences 
with lattice Monte Carlo data [TTl [T^] may seem more 
surprising, since in the particular case of the balanced 
system these algorithms are free of the sign problem, al- 
lowing in principle to approach the balanced unitary gas 
model in an unbiased way. However, eliminating system- 
atic errors from lattice-discretization and finite volume 
requires extrapolations which are either not done |llj 
or difficult to control [El [13]. The ENS experimental 
pressure EOS [55] lies systematically below ours, slightly 
outside the reported error bar. The experimental results 
from Tokyo [13] do not agree with the virial expansion at 
high temperature. The BDMC results agree excellently 
with the present experimental data all the way down to 
the critical temperature for superfluidity (see Fig. 4). On 
approach to (/3/x)c, we observe the growth of the corre- 
lation length in the BDMC pair correlation function F. 
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FIG. 4. Equation of state of the unitary Fermi gas in the normal phase, (a) Density n and (b) pressure P of 
a unitary Fermi gas, normalized by the density no and the pressure Po of a non-interacting Fermi gas, versus the ratio of 
chemical potential jj, to temparature T. Blue filled squares: Bold Diagrammatic Monte Carlo (this work), red filled circles: 
MIT experiment (this work) . The BDMC error bars are estimated upper bounds on systematic errors. The MIT error bars are 
one standard deviation systematic plus statistical errors, with the additional uncertainty from the Feshbach resonance position 
shown by the upper and lower margins in red solid lines. Black dashed line and red triangles: Theory and MIT experiment (this 
work) for the ideal Fermi gas, used to assess the experimental systematic error. Green solid line: third order virial expansion. 
Open squares: first order bold diagram jI5.ll9| . Green open circles: Auxiliary Field QMC [TT]. Star: superfiuid transition point 
from Determinental Diagrammatic Monte Carlo [13]. Filled diamonds: ENS experimental pressure EOS [22]. Open pentagons: 
pressure EOS from Tokyo experiment |23] . 



A protocol for extracting the critical temperature itself 
from the BDMC simulation will be presented elsewhere. 

We are not aware of any system of strongly correlated 
fermions in Nature where experimental and unbiased the- 
oretical results were compared at the same level of accu- 
racy. Even for bosons, the only analog is liquid ^He. This 
promotes the unitary gas to the major testing ground for 
unbiased theoretical treatments. The present BDMC im- 
plementation should remain applicable at finite polariza- 
tion and/or finite scattering length, opening the way to 
rich physics which was already addressed by cold atoms 
experiments [10, 24 27J. We also plan to extend BDMC 
to superfiuid phases by introducing anomalous propaga- 
tors. Moreover, since the method is generic, we expect 
numerous other important applications to long-standing 
problems across many fields. 

Note added in proof. After a preprint of this work be- 
came available, new Auxiliary Field QMC data were pre- 
sented [3T], with undetermined systematic errors whose 
evaluation in future work is called for by the authors of 
Ref. m]. 



Methods 

The experimental setup has been described previ- 
ously [23]. In short, ultracold fcrmionic ®Li is brought to 
degeneracy via sympathetic cooling with ^^Na. A two- 
state mixture of the two lowest hyperfine states of ^Li 
is further cooled in a hybrid magnetic and optical trap 
at the broad Feshbach resonance at 834 G. We employ 
high-resolution in situ absorption imaging to obtain the 
column density of the gas, that is converted into the full 
3D density via the inverse Abel transform |28^ . Equiden- 
sity lines provide equipotential lines that are precisely 
calibrated via the known axial, harmonic potential (axial 
frequency i^^ — 22.83 ± 0.05 Hz). Equipotential averag- 
ing yields low-noise profiles of density n versus potential 
V. Density is absolutely calibrated by imaging a highly 
degenerate, highly imbalanced Fermi mixture, and fit- 
ting the majority density profile to the ideal Fermi gas 
EOS ^4 . In contrast to previous studies [211 1M| , our 
analysis does not rely on the assumption of harmonic 
trapping. 

Thermometry is performed by fitting the density pro- 
file to the EOS constructed thus far, restricting the fit to 
the portion of the density profile where the EOS is valid. 
In the high-temperature regime, the EOS is given by the 
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virial expansion 

nX^ = + 26262'^'^ + 8636^^^ + ... (1) 

where the virial coefficients are &2 = 3\/2/8 [29], and 
&3 = -0.29095295 Fitting a high-temperature 

cloud to the virial expansion gives the temperature T 
and the chemical potential /i of the cloud, and the EOS 
nX^ — f{f3fi) can be constructed. We have used Eq.([l]) 
for /3/i < {f3fi)max — —1.25 and we checked that our EOS 
did not change within statistical noise if we instead used 
iPlJ')max = —0.85. Once a new patch of EOS has been 
produced, it can then in turn be used to fit colder clouds. 
Iteration of this method allows us to construct the EOS 
to arbitrarily low temperature. A total of ^ 1000 pro- 
files were used, with 10 to 100 profiles (50 on average) 
contributing at any given 
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